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ABSTRACT 

We analyze the time evolution of the luminosity of a cluster of Population III protostars 
formed in the early universe. We argue from the Jeans criterion that primordial gas can col¬ 
lapse to form a cluster of first stars that evolve relatively independently of one another (i.e., 
with negligible gravitational interaction). We model the collapse of individual protostellar 
clumps using 2+1D nonaxisymmetric numerical hydrodynamics simulations. Each collapse 
produces a protostar surrounded by a massive disk (i.e., M f ]j s k/A/, > 0.1), whose evolution 
we follow for a further 30-40 kyr. Gravitational instabilities result in the fragmentation and 
the formation of gravitationally bound clumps within the disk. The accretion of these frag¬ 
ments by the host protostar produces accretion and luminosity bursts on the order of 10 6 L 0 . 
Within the cluster, we show that a simultaneity of such events across several protostellar clus¬ 
ter members can elevate the cluster luminosity to 5-10x greater than expected, and that the 
cluster spends ~ 15% of its star-forming history at these levels. This enhanced luminosity 
effect is particularly enabled in clusters of modest size with ~ 10-20 members. In one such 
instance, we identify a confluence of burst events that raise the luminosity to nearly 1000 x 
greater than the cluster mean luminosity, resulting in L > 10 8 L 0 . This phenomenon arises 
solely through the gravitational-instability-driven episodic fragmentation and accretion that 
characterizes this early stage of protostellar evolution. 

Key words: accretion disks—cosmology: theory—hydrodynamics—stars: clusters—stars: 
formation—stars: Population III 


1 INTRODUCTION 


Following the emission of the Cosmic Microwave Background, 
the universe entered a cosmological “dark age” that ended with 
the formation of the first stars in the universe. These first stars 
(known as Population III stars) were responsible for producing 
the u ltraviolet radiation that be gan the reionization of the universe 
(e.g jTumlinson & Shullfcoooh . and their supernovae were respon¬ 
sible for enrichi ng the intergalactic medium with the first heavy 
elements (e.g., Miralda-Escude & Reesl 1 19971 : iGnedin & Ostrikea 
ll997l:lFerrara et al.ll2000h.~ 

Cosmological-scale simulations that follow both the dark mat¬ 
ter and baryonic components of the early universe have yielded the 
consensus opinion that Population III stars formed in dark matter 
minihalos with masses of approximately 10 lj M 0 . These 3 <r+ per¬ 
turbations over the backgro und dark matter dens i ty field virialized 
by redshifts o f z ~ 20-50 dTegmark et al.lll997l :jAbgletal.J|2002l: 
IBromm et al.ll2002l) . With few exceptions (e.g.. iTurk et alj|2009h . 
these simulations suggest that the gas pooling into the halos un¬ 
derwent a quasi-hydrostatic contraction until t hey had sufficient 


mass to trigger run a way gravitational colla pse dAbeLetafl 


Bromm_eUd. 200i_ Bromm & Loeb|_ 2004 lYoshida et al 


2002 ; 


2006; 


O’Shea & Normanl2007l : lYoshida et alj2008h . These studies estab 
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lished the standard paradigm that the progenitor cloud cores of the 
first stars were most likely to have been massive and formed in rel¬ 
ative isolation. 

In contrast, it is well understood (theoretically as well as ob- 
servationally) that most star formation in the present-day universe 
arises from the fragmentation of molecular clouds, resulting in a 
multiplicity of young st ellar objects being formed in close prox- 
imity_to each other (e.g., ICarnenter et alJfl997l : Iflillenbrandll 1 9971 : 
lLada & Ladall2003h . Motivated by this, several recent studies have 
explored the fragmentary nature of primordial gas in the early uni¬ 
verse, and have been able to resolve fragmentation in the disk-like 
environments surro unding the first protostars, thus cha l lenging the 
stand ard paradigm dStacy e t alj|20ld : Idarket all201ll : iGreif et afl 
l201ll : lvbrobvov et al.ll2013l . hereafter VDB 2013). 

Clearly some ambiguity remains regarding the initial condi¬ 
tions and the formation mechanism(s) of the first stars. Observa¬ 
tions will be required to accurately distinguish between the many 
existing theories of Population III star formation and evolution. In 
fact, the detection of primordial star clusters and galaxies in the 
early universe has already been defined as a majo r goal for next 
genera tion telescopes (e.g.. Iwindhorst et al.l 120061) . IBromm et al.l 
( 2001 1 were among the first to investigate the spectral energy dis¬ 
tribution of primordial stars theoretically. Later studies have ex¬ 
panded on their results to show that isolated Population III stars are 
likely to be too faint for detection by instruments such as the forth- 
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coming James Webb Space Telescope, even when their fluxes are 
enhanced via chance gravitational lensing dRvdberg et al.lljoTjl) . 
Several authors have also turned their attention toward the poten¬ 
tial of observing clusters, dwarf galaxies, and massive galaxies 
that contain Population III stars that may have formed at lower 
redshifts (i.e., 2 < 10) due to inhomogeneous metal enrichment 
of the intergalactic medium following the first supemovae (e.g.. 


Ciardi & Ferrara 200ll; Scannanieco et al. 

200ji Tornatore et al. 

20071; Johnson 2013; 

Safranek-Shrader et al 

2014K. 


In this paper we propose a scenario for the formation of a clus¬ 
ter of Population III stars. We argue that the gas pooling into the 
dark matter halos in which the first stars formed is subject to the 
Jeans criterion analogously to the fragmentation of giant molec¬ 
ular clouds in the present-day universe. As a result, these halos 
were capable of producing small clusters of first stars. Using non- 
axisymmetric numerical hydrodynamics simulations, we study the 
gravitational-instability-driven fragmentation and accretion in the 
collapsing protostellar environment. The resulting burst mode of 
accretion is even more prominent in a Population III environment 
than in present-day star formation, as shown by IVDB 2013 . We 
use the calculations for individual cluster members to compile the 
frequency, magnitude, and luminosity of burst events for each clus¬ 
ter as a whole. We find that a simultaneity of accretion events can 
produce bursts of luminosity that are several orders of magnitude 
greater than the mean cluster luminosity. 

The structure of this paper is as follows. In Section 2 we de¬ 
scribe how gas that settles into the host dark matter halos is subject 
to the Jeans fragmentation criterion, allowing for the formation of 
a cluster of first stars. In Section 3 we describe our numerical sim¬ 
ulations (as well as our selections for the initial conditions) for the 
formation of each protostar, and calculate the luminosity for each 
of these cluster members. In Section 4 we discuss the implications 
of having a multiplicity of protostars simultaneously experiencing 
bursts of accretion, and calculate the effect this has on the lumi¬ 
nosity of the cluster. We also discuss the implications of this phe¬ 
nomenon for future observational programs. Finally, in Section 5 
we conclude with a brief discussion of our results. 


2 THE CASE FOR POPULATION III STAR CLUSTERS 

The formation of a cluster of Population III stars is thought to 
begin with the collapse of low density baryon gas into the gravi¬ 
tational potential wells established by dark matter minihalos that 
have collapsed and virialized by z ~ 20-50. The gas to dark mat¬ 
ter fraction in these halos is roughly 10%, and amounts to a gas 
mass of 10 4 to 10 s Mg. Numerical simulations of this collapse re¬ 
veal that the gas streaming into th ese potential wells exhibits fila¬ 
mentary and knotty stru cture (e.g.. lBromm et ali i 1999: Cla rk et al.l 
l2008l:lGreifeTai]|201ll) . Indeed, gravity is well known to enha nce 
such anisotropic structure during collapse (e.g.. lLin~et all 19651) . 

The formation of H 2 cools the gas efficiently (and lowering 
the Jeans mass) from temperatures of a few times 1000 K to a few 
times 100 K, allowing the gas density to increase to n ~ 10 4 cm -3 . 
While the imprint of substructure exists within the gas morphology, 
the formation of H 2 is inefficient for cooling the gas below a tem¬ 
perature of a few times 100 K, inhibiting further collapse. Instead, 
these precursor imprints of fragmentation must next undergo a slow 
quasi-hydrostatic contraction as they a c crete additional mas s (e.g., 
iTegmark et alJll997l : lAbel et alj 120021 : iBromm et alj|2002h . Run¬ 
away gravitational collapse is only then triggered when the mass 
of these “weak clumps” exceeds the local Jeans value at this scale, 


being (e.g. Jdarke & Brommll2003l) 




. -1/2 


( 1 ) 


For example, a halo with a gas mass of 10 5 M 0 and a 10% star 
formation efficiency would form a weak cluster with ~ 20 mem¬ 
bers. Though this sequence of events differs slightly from that of 
present-day star formation, the latter also envisions a multiplicity 
of approximately Jeans mass fragments that form in close prox¬ 
imity to pr oduce a weak or strong cluster (e.g., in Taurus and p 
Ophiuchus; lOnishi et al.lll998l ; ljohnstone et al.ll2000l) . 

These resultant massive clumps, each containing roughly 
400 Mg of gas, are the sites of first star formation within the halo. 
The dynamical state of a typical clump formed in this way—its 
mass, physical extent, temperature, and angular momentum—are 
all determined by the collapse process. The initial conditions for 
the further contra ction of these clumps are t herefore relatively well 
constrained (e.g jYoshida et alj2003l 20061) . 

Authors such as Abel et all d 19981) have estimated that the ef¬ 
ficiency of this fragmentation—with which the gas pooling into the 
dark matter halo is assimilated into high-density clumps—could 
vary between as little as a few percent to nearly 50%. As the larger 
clumps tend to retain their individuality against dispersal and/or 
mergers, the final evolutionary state of the halo gas is expected to be 
a small cluster of isolate d objects (e.g. JClark et al.l2008l : lTurk et al.l 
120091: iGreif etai1l2012l) . However, owing to lingering ambiguities 
about the manner in which this occurs, we adopt the definition of 
a “cluster” as being any association of 2+ stars forming from in¬ 
dependent gaseous clumps that result from the subfragmentation 
of primordially pristine gas that has pooled into a single dark mat¬ 
ter halo. Additionally, we concern ourselves with only the forma¬ 
tion period during which the cluster members have masses below 
40 Mq. This allows us to assume that the predominant gas fraction 
within the halo is not ionized by the star formation as the cluster 
evolution proceeds. In fact, such conditions represent the analogue 
of present-day so-called “embedded” clusters, in which more than 
80% of the cluster members belong to the Class II/III evolution¬ 
ary phases, and the mass function of the cl uster is assumed to be 
no longer evolving (e.g. JGutermuth et alj|2009h . We also assume 
that each clump is able to promptly form stars within ~lMyr, 
or roughly the lifetime of the most massive individual stars (e.g., 
iBond et all 1981) . 

Some caveats to our model assumptions are that we are as¬ 
suming cluster formation according to the standard hydrodynamic 
Jeans criterion, and disk evolution in the hydrodynamic limit. Clus¬ 
ter formation may be strongly affected by several other factors. 
Turbulence generated by colliding flo ws may leave a strong initia l 
imprint that dominates fragmentation dHeitsch et alJl2006l . 120081) . 
Radiative feedback from the very first massive stars can also af- 
fect the further fragmentation of a cloud or disk ( IStacv et al.l2012l) . 
Furthermore, if a seed magnetic field exists and is amplified by dy¬ 
namo action during primoridal cloud collapse dTurk et al.ll2012l) . 
then magnetic field effects can signific antly affect the Jeans scale 
dCiolek & B as ul 20061: BasuetaL 2009) and affect global mass flow 
in the cloud dPrice & Batel 120081) . Magnetic fields can also in¬ 
hibit the formation of large disks or their fragmentation in sim 


ulations of present-day star formation (e.g., [M ellon & Li 
Commen^on'etalj^ 20 id : iMachida et alJ 1201 if : IPanp et all 


2008 


2012 


Tomida et alj 2015 ). and their effect in the Population III regime 


remains to be further elucidated. 
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3 INDIVIDUAL CLUSTER MEMBERS 

The clumps that emerge from the fragmentation are assumed to be 
relatively isolated. Competition for accretion between clumps (due 
to protostellar crowding) and effect s arising from gravitational in¬ 
teract ions are negligibly small (e.g. jTurk et al]l2009t iHirano et al.l 
|20I4|) . We thus model the formation and evolution of each clus¬ 
ter member with its own unique simulation—the initial conditions 
of which stem from the arguments of the preceding section. Here 
we provide details about the the numerical aspects of our code, our 
specific choices of initial conditions, and investigate the long-term 
behavior exhibited in a fiducially constructed model. 


3.1 Numerical Simulations 

We carry out numerical simulations of the gravitationally in¬ 
duced collapse of the primordial gas in 2+ID, assuming a thin- 
di sk geometry. Our code is a mo dified version of that presented 
in IVoro bvov & Basu (2005. 12 OO 6 I) . The hydrodynamic equations 
are solved using a finite differ ence scheme with a tim e-explicit 
operator-split solution, based on lstone & Norman! il992f t. A thor¬ 
ough description of our code is presented in lVDB 20131 . 

The mass and momentum transport equations in the thin disk 
limit can be expressed as 

r)Y 

— + V-(Ev) = 0, (2) 

(Ev) + V ■ (E v O v) = -VP + E^, (3) 

ot 

in which E is the surface mass density, v is the velocity of the disk 
material, P is the vertically integrated pressure (determined assum¬ 
ing the disk is in vertical hydrostatic equilibrium at all times), and 
V is the planar gradient operator. The gravitational acceleration in 
the plane of the disk ( g) includes the contribution from the protostar 
(once formed), from material inside the sink cell, and from the self¬ 
gravity of the disk and surrounding cloud core. All vectoral terms 
and quantities are understood as having only r and cf> components 
in this formulation. 

In order that the environment of the collapsing core accurately 
reflect primordial conditions, equations 0 and l[3j are closed with 
a barotropic_relation that fits the ID core collapse simulations of 
[Omukai et al.l (l2005h . These include the detailed chemical and ther¬ 
mal processes of the collapsing gas. Additional details of our model 
are provided in lVDB 200 . 


3.2 Initial Conditions 


The initial conditions for the radial gas surface density S and 
angular velocity 9 profiles for primordial cores are taken to be 
very similar to t hose o f present-day star formi n g cores (e.g., 
Omukai & Nishil 19981: Vorobyov & Basul 120061 : lYoshida et al.l 


20081: 1 Vorobyov & Basdl2010b. 


S = E 0 



n = 20o 



(4) 

(5) 


The radial surface mass de nsity profile E is that of an integrated 
Bonnor-Ebert sphere jDapp & Basul 120090 . while the form of the 


initial angular velocity profile H corresponds to the differential ro¬ 
tation profile expected for a cor e c ollap sing out of a near-uniform 
initial surface density field dBasull 1997h . 

We constrain Eo by assuming a constant ratio of the radius 
of the outer computational boundary r out to that of the centrally 
plateaued region ro: r ou t/ro = 6; so that each core has a simi¬ 
lar initial form with Eo « 0.25 gcm~ 2 . The parameter r ou t thus 
also determines the mass of the core, which for r ou t = 0.5 pc is 
approximately 300 Mq. These choices are constant for each of the 
cores simulated, and is consistent with the typical size and mass 
found from ab initio cosmological simulations of the collapse of 
primordial starless cores (e.g., Yoshida et alhoQfj ). 

The angular momentum of the cloud core is parameterized 
by Ho—as appears in equation 0 — and is related to the dimen¬ 
sionless parameter 7 ] = (Horo/cs) 2 teas 3 11997b . ri is related 
to the ratio of the rotational to gravitational energy of the cloud 
core, /3 = E ro t/\E g \, with /3 = 0.9r;. Finally, /3 relates to the 
so-called “spin parameter,” A = \ffi, that is most often used 
to characterize the angular momentum of dark matter halos (and 
their associated gas) formed from cosmologica l initial conditions 
(e.g., iBames & Efstathioulll987l : lRvdenlfl988b . The spin param¬ 
eters of each cloud core are lognormally distribut ed w ith mean 
A = 0.05 and variance er 2 = 0.5 (following IClardnetl 1200 ll : 
lo’Shea & Normanil2007h . 

Each model is run on a polar coordinate grid with 512 x 512 
spatial grid zones in r and cj>. The inner and outer boundary con¬ 
ditions allow for free outflow from the computational domain. Ra¬ 
dial grid points are logarithmically distributed to allow for better 
numerical resolution toward the innermost region of the disk: the 
innermost cell outside of the sink region has a radius of approxi¬ 
mately 0.1 AU and is 1.9 AU (both radially and azimuthally) at a 
radiu s of 10 0 AU. 

iTan & McKeel d2004b and lHosokawa et all <1201 ll) have shown 
that beyond 30-40 Mq the role of increasing stellar luminosity be¬ 
comes critically important to understanding the subsequent evolu¬ 
tion of these protostars, as the intensely ionizing radiation begins 
to inhibit H 2 formation, which is the primary coolant in the gas. 
We therefore terminate our simulations once the protostar reaches 
amass of 40Mg. 

3.3 Evolution of the Protostellar Disk 

Our fiducial model is characterized as a clump of gas roughly 0.5 pc 
in radius, with a mass of ~300 Mq, and at a temperature of 300 K. 
The spin parameter A of the clump is equal to the mean of the dis¬ 
tribution, A ~ 0.05. The collapsing cores we model compare well 
to those de rived from 3D numerical simulations inside primordial 
minihalos dYoshida et alj|2003 : Idark et al.ll201 lb . Differences are 
attributable to our cores being toward the lower end of the mass 
spectrum that has been studied by these authors. 

A quasi-Keplerian disk forms around the protostar within 
~3kyr of the formation of the central protostar (in our fiducial 
model). The disk begins to fragment within a few hundred years af¬ 
ter itsformation/This time scale is somewhat longer than that found 
bv Idark et al.l d201 lb and Icireif et aTI d20 1 ll ). who used sink cells 
with smaller radii of 1.5 AU. However, this timescale is similar to 
that found by lSmith et al.1 d2012ab . who used sink cells with com¬ 
parable radii of 20 AU (the central sink cell in our simulation being 
~10 AU in radius). 

In Figure|T]we present three snapshots of the disk surface den¬ 
sity inside of a 1000 AU radius and spanning lOkyr of the disk 
evolution. The left panel shows that a rich density structure already 
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Figure 1 . Projection of the disk surface mass density S within a 2000 X 2000 AU volume centered on the accreting protostar. The time in each frame is t = 5, 
10, and 15 kyr after formation of the protostar. Left: Some regions of the disk are Toomre-unstable, and several fragments can be identified at radii throughout 
the disk, between 10 and several hundred AU. Center: As these clumps are accreted onto the protostar or dispersed, the disk enters into a moderately quiescent 
period during which it is relatively stable against fragmentation. No fragments are present in the inner 300 AU of the disk. The three large fragments at radii 
beyond this are remnants from an earlier phase of fragmentation that had been raised to higher orbits. Right: Continued accretion from the parent cloud core 
builds up the mass of the disk, eventually triggering another period of vigorous fragmentation. The two fragments immediately to the left of the sink cell in this 
panel provide an example of how larger fragments can be sheared apart prior to being accreted through the sink cell, leading to intense and rapid variability in 
the accretion luminosity of the protostar. 


exists within the disk a mere 5 kyr after formation of the central 
protostar. Together with the right panel, the clearly defined con¬ 
densations of gas that have fragmented out of the disk indicate that 
the disk experiences multiple episodes of instability. The central 
panel presents an intermediary quiescent phase during which there 
is little to no fragmentation and the accretion of gas onto the pro¬ 
tostar occurs relatively smoothly. The number of fragments within 
the disk clearly varies with time, as some fragments are tidally dis¬ 
persed and others are accreted onto the central protostar; a result of 
the gravitational torques exerted by the fragments on each other and 
from larger spiral arm structures that form within the disk. Though 
accretion gradually drains the disk, new episodes of fragmentation 
are continually stimulated by the resupply of fresh gaseous ma¬ 
terial from the surrounding core envelope. It is also evident that 
most of the fragments are formed in the intermediate to outer disk 
region (> 50 AU), which is consistent with the numerical simula¬ 
tions of disks around protostars in the present-day universe (e.g., 
IStamatellos & Whitworthll2008l : IClarkell2009l) . 

Those fragments that pass through the inner computational 
boundary are assumed to be readily accreted by the central proto¬ 
star. Though this boundary is not the protostellar surface itself, we 
can use the temperature (or more appropriately, the sound speed) of 
the disk material along this boundary to estimate the instantaneous 
mass accretion rate onto the protostar as 


M : 


Cs_ 

G 


( 6 ) 


( Shull 197 7!). In our fiducial model we find c s is about 2.0kms _1 . 
This corresponds to an instantaneous mass accretion rate of 
~10 ~ 3 M@ yr™ 1 . We note that this is just below the critical (Ed¬ 
dington) accretion rate of MecM — 4.0 x lO™ 3 M 0 yr _1 above 
which the corresponding radiatio n pressure alone is capable o f halt¬ 
ing the accretion flow entirely (Hoso kawa & Qmukail l2009). This 
estimate matches the time averaged value of the mass accretion 
rate experienced by our model protostar (Figure |2), and is consis¬ 
tent with several other models of the runaway collapse of primor¬ 
dial gas as has been determined both analytically and from simu- 
lations ( Omukai_&Palla[ 2003; Bromm&Loeb| 120041 : Idark et ail 
1201 ll:lHosokawa et al.11201 1 ; Smith et al.ll2012bl). 


We present the complete mass accretion history for our fidu¬ 
cial model in the top panel of Figure [2] The formation of the pro¬ 
tostar is marked by the sharp rise in the mass accretion rate, which 
we define as time t — 0. Following this, material from the sur¬ 
rounding envelope of the progenitor cloud core continues to stream 
onto the protostar at a rate of a few times 10™ 3 Mg yr -1 . The mass 
of the protostar rapidly increases (within ~4 kyr) to approximately 
15 Mg before the accretion flow is temporarily halted by the for¬ 
mation of a quasi-Keplerian disk. 

The disk modulates the subsequent accretion of disk mate¬ 
rial onto the protostar, with gravitational torques redistributing the 
mass_and angular momentum of the infalling cloud core material. 
In VDB 2013 we demonstrated that the disk self-gravity quickly 
induces the formation of spiral arms. Furthermore, sections of the 
disk in which the local value of Toomre’s Q falls below unity be¬ 
come subject to fragmentation dToomret 19641 ). These fragments are 
then torqued inward along ballistic trajectories before ultimately 
being accreted by the protostar, resulting in the episodic bursts of 
accretion (see Figure[2]i. 

The cumulative effect of a quiescent mode of accretion that is 
punctuated by the episodic bursts is also evident in the curves of 
mass growth in the bottom panel of Figure[2] The protostellar mass 
(shown in blue) grows rapidly during the initial phase of smooth 
accretion (which lasts ~4 kyr). However, the burst mode of accre¬ 
tion comes to dominate the subsequent growth of the protostar as 
is evident by the abrupt increases in the protostellar mass, typi¬ 
cally a few Mg at a time. These increases are typically followed by 
plateaus during which the mass of the protostar changes very little 
as the disk equilibrates back into a quasi-Keplerian state. Each burst 
event is mirrored by a corresponding decrease in the total disk mass 
(in orange). However, the overall mass of the disk actually contin¬ 
ues to increase in time due to the steady accretion of material from 
the remnant of the progenitor cloud core (the dashed black line in 
the top panel of Figure |2j. 

Additional details of the effects of the gravitational torques on 
the disk, and the resulting recurrent character of the disk fragmen¬ 
tation, are discussed in VDB 2013. Here we focus on the signature 
of this behavior on the protostellar accretion luminosity. 
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Figure 2. Top: Temporal evolution of mass accretion rates: from the cloud 
core onto the disk at 3000 AU (dashed black), and from the disk onto the 
protostar (solid black). The rectangular window time-averaged (in 1 kyr 
intervals) mass accretion rate—the quiescent accretion rate—is in orange. 
Bottom: Growth of the protostellar (blue) and disk (orange) masses in time. 
The protostellar mass grows via punctuated equilibrium, while the episodic 
increases in the mass of the protostar are coincident with the episodes of 
decline in disk mass. Nevertheless, the disk mass increases with time owing 
to the continually replenishment of disk material from the surround core 
material. 


3.4 Accretion Luminosity 


A protostar’s luminosity is a product of competition between mass 
growth from accretion and radiative loss from the protostellar inte¬ 
rior. However, it is not until the protostar begins contracting toward 
the main sequence that its internally generated luminosity L* sur¬ 
passes the accretion luminosity L acc • During the earliest stages of 
its evolution, the source of a protostar’s luminosity is almost en¬ 
tirely fr om accretion, so that L ~ L acc (up to masse s o f M* ~ 30- 
40 Mr. jTan & McKedbOO^fHosokawa et all201lh . In lVDB 20131 
we showed that the large variability in accretion experienced by the 
protostar results in accretion luminosities several orders jrf magni- 
tude greater than might other wise be expected (comp are ClarketaU 
1201 ltlHosokawa et alj|201 ll. and lSmith et ai]l2012bl to lVDB 20131 


for example). 

We estimate the accretion luminosity assuming that any ma¬ 
terial landing on the surface of the protostar has its kinetic energy 
dissipated radiatively at a rate 


r _ GM t M 

L -^mT 


(7) 


where R * is the protostellar radius. 

In the absence of a detailed model for the stellar interior we 
instead fit the evolutionary models of lOmukai & Fallal d2003i) with 
a piecewise power-law approximation to the radial expansion of the 
protostellar surface as a function of the protostar mass. Following 
formation of the hydrostatic core, the protostellar radius is expected 
to grow according to a mixed power-law as R* ocM *' 27 M^ 1 ; 
where M_ 3 denotes the ratio of the actual instantaneous m ass ac¬ 
cretion rate M t o a va lue of 10 -3 M@ yr -1 JStahler et alJll98fi| : 
lOmukai &~P alla 2003). Increasing temperature within the core 
drives a luminosity wave outward, causing a rapid expansion of 
the stellar surface. When this wave breaches the protostellar sur¬ 
face, the interior is able to relax and the protostar begins Kelvin- 
Helmholtz contraction toward the main-sequence. 

The following relations approximate the evolution^ of R » 
through its transitions through these phases dSmith et~al]|20 1 2rj) : 

r 26 M*<M 1 

R ., = ( Mf Mi < M, < M 2 ■ (8) 

( A 2 MP 2 , M* ^ M 2 and R* < Rms 


The constants Ax and A 2 are matching conditions that ensure the 
functional form of R * is smoothly varying between transitions. 
The mass parameter Mi marks the transition between the adia¬ 
batic phase of growth and the arrival of the luminosity wave at 
the protostellar surface; M 2 , the transition between the luminosity 
wave driven expansion and subsequent Kelvin-Helmholtz contrac¬ 
tion. Mi and M 2 are fixed by the instantaneous mass accretion rate 
as the protostar transitions between phases, and are defined as 


Mi = 5 M°-f 7 , 
M 2 = 7M°3 7 . 


(9) 


Note that as the evolution of the protostellar interior occurs roughly 
adiabatically, due to the long cooling time therein, the variability in 
the accretion rate induced by the burst mode of accretion does not 
result in significant variability in the radius of the protostar. 

The accretion luminosity calculated for our fiducial model is 
shown in Figure[3](black). In overlay is the time averaged rate (i.e., 
the quiescent rate, L q , in orange), and a demarcation 10 x greater 
than the quiescent rate (dashed orange). Accretion events during 
which L exceeds L q by the factor of 10 are highlighted (dashed 
blue). Once the protostar is formed (designated as t — 0), L climbs 
very quickly to ~10 4 Lq. The luminosity remains at about this 
level during the period of smooth accretion until a disk is formed 
at £ ~ 4 kyr. Although the mean rate remains at roughly 10 4 Lq, 
the episodic nature of the subsequent vigorous accretion gives rise 
to significant variability, with some peaks reaching several times 
10 6 Lq —two orders of magnitude greater than L q . One might ex¬ 
pect this large variability in luminosity caused by the individual 
burst events to affect the accretion flow via feedback from the en¬ 
hanced radiation field. However, as the duration of the individual 
bursts are quite short (typically on the order of a few times 100 yr), 
we expect that they do not have any appreciable effect on the long¬ 
term growth of the protostar. 

In the bottom panel of Figure [3] we focus on a 1 kyr win¬ 
dow of our simulation in order to highlight our burst identification 
scheme. We calculate the effective mean accretion rate (the solid 
orange line) using a moving 1 kyr window. As burst events are ul¬ 
timately attributable to the accretion of large individual fragments 
(with masses on the order of 0.1-1 Mg) and have an overall mean 
duration of a few times 100 yr, the use of a 1 kyr moving window 
allows us to supress variations due to the large fragments being 
sheared apart into several smaller ones, each of which is accreted 
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Figure 3. Top: Temporal evolution of the (accretion) luminosity from the 
protostar in our fiducial model (in black). In overlay is the time averaged 
quiescent rate (solid orange), and the 10xL q rate (dashed orange), with 
those accretion events whose peak luminosity exceeds this threshold high¬ 
lighted in blue. Bottom: A zoom in on a 1 kyr window of the luminosity 
history for clarity. Accretion events in excess of 10xL q are clearly identi¬ 
fied (labeled, a, b, and c, and highlighted in blue). 

by the protostar in rapid succession ( VDB 20 I 3 I) . Hence, we also 
include as part of a single burst event all points left- and rightward 
of the peak luminosity that are greater than the mean (and not only 
those points for which L is strictly > 10xL q ). 


4 CLUSTERS OF POPULATION III PROTOSTARS 

Figure[4]illustrates how the burst mode in individual protostars can 
sum to produce a cluster luminosity. For simplicity, we consider a 
cluster of only two stars, with one being our fiducial model with 
A = 0.05, and the second characterized by A = 0.065 but forming 

5 kyr after the first. The panels in the top- and middle-left column 
are the luminosities of the individual cluster members while the 
bottom-left shows the cumulative luminosity for the cluster as a 
whole. The histograms on the right-hand side present statistics of 
the fractional number distribution of burst events /b of different 
duration Afb- These can be used as metrics to evaluate differences 
between the individual and cluster luminosities. 

Prior to t = 5 kyr, before the second cluster member forms, 
the luminosity of the cluster is of course just the luminosity of the 


first and as yet sole member. The second cluster member begins 
its evolution from t ~ 5 kyr (where time t = 0 marks the for¬ 
mation of the first cluster member). In the subsequent ~2kyr, as 
the second cluster member smoothly accretes material from its sur¬ 
roundings, its contribution to the total cluster luminosity increases. 
In fact, comparing the cluster luminosity to the individual compo¬ 
nent luminosities reveals that a significant amount of the variability 
in the luminosity of the first cluster member is being completely 
obscured. During this time, the cluster luminosity is consistently 
above ~10 4 L©. By t ~ 7 kyr, both cluster members harbor their 
own protostellar disks, and accretion onto each protostar is being 
driven by the action of gravitational torques in the massive disks 
surrounding each host—i.e., via the burst mode of accretion. Cor¬ 
respondingly, the cluster luminosity thereafter exhibits a significant 
amount of variability, comparable to that of either of its component 
members. 

In the right-hand column of Figure [4] we present histograms 
of the duration of the burst events exhibited by the individual clus¬ 
ter members (top- and middle-right), and for the cluster as a whole 
(bottom-right). The first cluster member exhibits a total of 73 indi¬ 
vidual burst events over the course of ~20.5 kyr—the span of time 
from when the disk forms to when the protostar reaches a mass of 
~40 Mq and the simulation is terminated. The mean burst duration 
is found to be ~ 96.3 yr. This amounts to approximately 7.0 kyr 
that are spent exclusively in the burst mode of accretion, with the 
remaining 13.5 kyr spent in quiescent phases. By comparison, dur¬ 
ing the approximately 18.5 kyr from the time the second cluster 
member forms to when the simulations end, the second cluster 
member exhibits a total of 68 individual burst events. The typical 
burst duration is found to be ~80.2 yr, for a total of ~5.5 kyr spent 
accreting via the burst mode, and the remaining roughly 13.0 kyr 
spent in quiescent phases. 

From this we can conclude that the burst mode of accretion 
plays a significant role in the mass growth of Population III pro¬ 
tostars. In both cases about 30% of each protostar’s individual 
accretion history is spent in the burst mode—that is to say that 
gravitational-instability-driven fragmentary accretion is responsi¬ 
ble for one-third to one-half of the total accretion onto the first 
stars. This is actually a significantly higher proportion of time than 
has been observed in analogous simulations of present-day star for¬ 
mation, wherein the burst mode of accretion is thought to be re- 
sponsible for only ~10% of a protostar’s accretion history (e.g., 
IVorobvov & BasiJl200dl20l(il:IVDB 20131) . 

The number and frequency of the observed bursts is also af¬ 
fected by whether one considers each cluster member individually 
or simply considers the cumulative luminosity of the cluster as a 
whole. The first and second cluster member yield 73 and 68 burst 
events over 18.5 and 13.0 kyr, respectively—1 event roughly every 
200 yr in both cases. Performing the same analysis of the luminos¬ 
ity, but of the cluster as a whole, the perceived number of burst 
events would be only 74: far fewer than might be expected given 
the number of actual burst events experienced by the cluster’s com¬ 
ponent members. 

Hence, the quiescent mode is clearly the dominant mode of 
accretion. As it is unlikely for all members of a cluster to be dim¬ 
ming simultaneously, the variations below the quiescent rate in the 
luminosity of a single cluster member rate are largely obscured. 
The cluster luminosity is thus somewhat greater than might be 
suggested by a simple summation of the individual luminosities. 
Hence, the number of burst events as determined by examining the 
cluster luminosity is actually fewer than might be expected, be¬ 
cause individual burst events can only be identified when the lu- 
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Figure 4. Accretion luminosity L and distribution of burst durations of the individual cluster members (top and middle panels) and for their combined output 
(bottom panel) in our example N = 2 member cluster. The top panel is our fiducial model with A = 0.050; the middle panel is a model with A = 0.065. 
Left: Accretion luminosities L are in black; the quiescent luminosity L q is in orange; the dashed blue line indicates a luminosity 10xL q . Note that the top 
panel is the fiducial model also shown in Figure [3] Right: Fractional number distributions (/b) of burst durations (Atb) f° r the individual cluster members 
(top and middle), and for the cumulative cluster profile (bottom). The total number of identified bursts (JVjb) is indicated in the upper-right of each panel. The 
approximate mean burst duration observed for the individual cluster members are 96.3 and 80.2 yr, respectively. The mean burst duration for the cumulative 
cluster luminosity profile is ~72.0yr. 
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Figure 5. Normalized accretion luminosities 71/71 q (left column) and histograms of the fractional number distributions of burst durations (right column) for 
two additional clusters containing TV = 16 (top), and TV = 128 (bottom) members. The dashed blue lines in the left panels denote thresholds of 10 and 
100 X the quiescent luminosity; all other lines and colors are as they appeal' in Figure [4] The normalized quiescent luminosity in each case corresponds to 
approximately 4 X 10 5 and 4 X 10 6 Lq, respectively. During the evolution of the TV = 16 cluster a confluence of burst events at f, re 18kyr produce a 
particularly prominent luminous event that is 647 X more luminous than the cluster’s quiescent level. However, the competition between the increasing cluster 
quiescent luminosity L q and the potential for such overlapping bursts is evident in the reduced magnitude of the fluctuations above L q with increasing cluster 
size TV. In the right panels, the total number of identified bursts (TV^) is indicated in the upper-right. 


minosity of an individual star exceeds the luminosity of the cluster 
as a whole. As a result, the perceived duration of burst events also 
decreases. 

To compare the variations in the luminosities of clusters of dif¬ 
ferent size TV, we normalize each cluster luminosity L by its own 
quiescent rate L q , and in Figure [5] compare them for clusters with 
TV = 16 and TV = 128. Luminosities are plotted in black; quiescent 
cluster luminosities in orange; and the dashed blue lines denote lu¬ 
minosity levels 10 x and 100 X the quiescent rate L q . The panels 
show that for larger TV it is less likely that there are burst events 
for which the cluster luminosity reaches above a certain multiple 
of L q . Conversely, as TV increases, there is also the possibility that 
two or more cluster members will experience simultaneous inde¬ 
pendent burst events. We use Figure[6]to illustrate how the the time- 
averaged quiescent luminosity (L q ) increases with the number TV 
of cluster members. The best-fit line represents a linear relation¬ 
ship. For completeness, we model clusters as large as TV = 1000, 
however most estimates of Population III clustering would imply 
much smaller clusters, as discussed earlier in this paper. 

In Figure [ 7 ] the vertical axis shows the duration of time At 


as a fraction of the cluster’s total star-forming lifetime t that a 
cluster spends above a certain luminosity threshold. The horizon¬ 
tal axis shows the sizes TV G {2,4, 8, 16, 32, 64,128, 256,1024} 
of the individual primordial star-forming clusters. Each curve rep¬ 
resents a specific luminosity threshold: from top to bottom, of 
5,10, 20, 50, andlOO times each specific cluster’s quiescent lumi¬ 
nosity. The diamond symbols along a single value of TV thus indi¬ 
cate how much time that cluster spends at an elevated luminosity 
L/L q . As the typical time frame for an individual burst event is 
on the order of a few hundred yr, the likelihood of two or more 
protostars within a cluster expressing simultaneous burst events is 
quite low. As the number of cluster members increases however, 
such coincidences become increasingly likely. However, this effect 
is counteracted by the fact that each additional protostar in a clus¬ 
ter also contributes to increasing the mean cluster luminosity. This 
prohibits all but the most intense burst accretion events from being 
visible. A balance between these effects is achieved for relatively 
small clusters, TV ~ 16. For these cases, we find nearly 15% of 
the cluster’s star forming period is spent at a luminosity level 10 x 
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Figure 6. The time-averaged quiescent luminosity (L q ) in clusters of size 
N £ {2, 4, 8,16, 32, 64,128, 256,1024} (from left to right) represented 
by diamond symbols. (L q ) increases roughly linearly as a function of the 
number N of cluster members as shown by the best fit line. 



Figure 7. The fractional duration of total time, Af/f, that a a cluster of 
a given size N spends at an elevated luminosity. Shown are the values of 
Af/f for different values of elevated luminosity L/L q = 5,10, 20, 50, and 
100 (appearing sequentially top to bottom). Most curves shows a maximum 
at N ~ 16. Larger clusters exhibit decreasing variability above the mean 
quiescent rate ( L q ). 

greater than the mean rate, L q . This corresponds to an absolute lu¬ 
minosity approaching 10 s Lq. 


5 DISCUSSION & CONCLUSIONS 

We have presented a scenario for the assembly of a cluster of first 
stars formed from the pristine gas that is pooling into newly viri- 
alized dark matter halos at redshift 2 ~ 20-50. The Jeans cri¬ 
terion provides a context for the combination of self-gravity and 
anisotropic structure to give rise to fragmentation of the gas. We 
estimate that a typical dark matter halo with mass between 10 5 and 
10 6 Mq contains 10% by mass of gas, of which roughly 10% actu¬ 
ally goes into stars. The Jeans criterion then implies the formation 
of clusters containing tens of members (perhaps 10-50 protostars), 
though the specific numbers are difficult to estimate due to the am¬ 


biguity with which th e disp ersal and merger of clumps prior to their 
collapse occurs (e.g jGreif et alllOl 1 1 ). 

We employ numerical hydrodynamics simulations in the thin- 
disk limit to then investigate the long-term evolution of the indi¬ 
vidual clumps. We model their collapse self-consistently into the 
formation of the protostar and its surrounding disk. In each case 
a disk forms relatively quickly, regardless of the initial conditions, 
within a few kyr of the formation of the protostar. The disk is cen- 
trifugally supported, but its mass lends it to being gravitational un¬ 
stable. We observe a rapid and episodic fragmentation of the disk 
during which fragments having between 0.1-1.0 Mq are formed at 
typical radii of ~50 AU. The fragments are torqued inward toward 
the central protostar where they are then accreted, resulting in mas¬ 
sive bursts of luminosity that can exceed the mean rate by as much 
as two orders of magnitude, occasionally exceeding 10 6 Lq (as in 
Figure |4}. 

In the context of this formation scenario we analyze the lu¬ 
minosity profiles that are produced by a multiplicity of first stars 
having formed approximately coevally within a young cluster. We 
assume that each cluster member forms independently and harbors 
its own disk that is subject to gravitationally induced episodic frag¬ 
mentation and accretion (i.e., the burst mode of accretion). With 
increasing numbers of cluster members, the quiescent luminosity 
of the cluster is seen to steadily increase. However, increases in 
the number of cluster members also increase the probability that 
two or more members experience a burst of accretion simultane¬ 
ously. In one cluster simulation we observe one such particularly 
luminous super-posed accretion event during which the cluster lu¬ 
minosity rises to nearly 1000X it’s quiescent rate, to 2.48x 10 8 Lq 

(Figure[5)- 

Competition between these two effects results in clusters of 
a certain size {N ~ 16) spending a sizable fraction of their 
star forming life-time (roughly 15%) at luminosities lOx greater 
than might otherwise be expected with smooth accretion (e.g.. 
Smith et al.l2012al) . Each additional object added to a cluster raises 
the mean/quiescent luminosity of the cluster, and this in turn tends 
to suppress the prominence of individual and multiple burst events 
in clusters in which N > 16. Such moderately sized clusters are 
also the most easily manufactured, capable of being produced as¬ 
suming even a relatively low star formation efficiency, and are of a 
number as has been counted in a vari ety of theoretical Population 
III star forming scenario s to date (e.g. JStacv et alj20ld : IClark et al.1 
l201ll:lGreifetai1l201ll) . 

The possible formation of clusters of Population III stars in 
the dark matter halos of the earl y unive rs e prov ides a unique op¬ 
portunity for observations. iRvdberg et alj d2013l) estimated the to¬ 
tal luminosity from even the most upper mass estimates of Popu¬ 
lation III stars (i.e., those with masses ~300Mq) and found that 
even these were likely too faint to be observed by next-generation 
telescopes such as the James Webb Space Telescope. The likeli¬ 
hood for such instruments to observe even lower mass Population 
III stars is thus highly improbable. However, as we have shown, 
clusters of even low mass Population III stars (as others have also 
indicated may exist: Bromm et al. 2003; Clark et al. 2008; and Greif 
et al. 2012), are capable of producing luminos ities in e xcess of lone 
very massive Population III stars (for example fBromm et al.l <[2001 ) 
estimate luminosities in the range ~ 10 6 -10 7 Lq for masses in the 
range 100-500 Mq). Observational constraints on the early stages 
of reionization, from 21 cm observations by the Square Kilometer 
Array (e.g., Carilli et al. 2004) may actually provide indirect con¬ 
straints on the abundance of such clusters. Observations by the At¬ 
acama Large Millimeter/submillimeter Array (ALMA) of emission 
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from dusty galaxies in the early universe may also provide insight 
into their possible formation as influenced by the formation of an 
initial cluster of Population III stars. However, it will be in the next 
decade that next-generation telescopes may actually provide direct 
constraints on the Population III star formation taking place in the 
early universe; with clusters of Population III stars being some of 
the most luminous at that epoch. 
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